Introdução

“Infecção hospitalar é uma infecção adquirida após a admissão do paciente na unidade hospitalar e pode se manifestar durante a internação ou após a alta. Pela sua gravidade e aumento do tempo de internação do paciente, é causa importante de morbidade e mortalidade, caracterizando-se como problema de saúde pública.”

Os dados disponíveis no arquivo data/SCENIC.txt foram coletados pelo CDC (US Center for Disease Control), no âmbito do Projeto SCENIC. O principal objetivo do projeto era determinar se programas de vigilância e controle foram capazes de reduzir as taxas de infecção hospitalar. Os dados referem-se a uma amostra de 113 hospitais selecionados a partir de um conjunto de 338 hospitais avaliados. Cada linha do conjunto de dados contém uma identificação (1-113) e fornece informação a respeito de 11 variáveis para um único hospital. Os dados apresentados referem-se ao período de 1975-1976.

As 12 variáveis são:

  1. IDnumber: 1-113 (identificação do hospital)
  2. LengthStay: período de internação médio de todos os pacientes no hospital (em dias)
  3. Age: idade média dos pacientes (anos)
  4. InfectRisk: risco de infecção, calculado como a probabilidade média estimada de contrair infecção no hospital (em %)
  5. CultRatio: razão do número de culturas realizadas pelo número de pacientes sem sintomas de infecção hospitalar, vezes 100
  6. XrayRatio: razão do número de raios-X realizados pelo número de pacientes sem sintomas de pneumonia, vezes 100
  7. NBeds: número médio de leitos do hospital, durante o período avaliado
  8. MedSchool: afiliação a alguma Escola de Medicina (1=Sim, 2=Não)
  9. Region: região geográfica (1=NE, 2=NC, 3=S, 4=W)
  10. DailyCensus: número médio de pacientes no hospital por dia, durante o período avaliado
  11. NNurses: número médio de enfermeiros no hospital
  12. Facilities: percentual de 35 serviços providos pelo hospital

Acredita-se que o período médio de internação de um paciente LengthStay (variável de resposta) possa ser previsto a partir do risco de infecção hospitalar, bem como outras características do hospital e de procedimentos de rotina realizados.

Análise Exploratória de Dados

Conduza a análise exploratória da massa de dados SCENIC, a fim de compreender suas características principais.

rm(list=ls())
#Carregando os dados em uma tabela
scenic <- read.table("data/SCENIC.txt", header = FALSE)

#Verificando a estrutura dos dados e atribuindo os nomes adequados às variáveis
str(scenic)
## 'data.frame':    113 obs. of  12 variables:
##  $ V1 : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ V2 : num  7.13 8.82 8.34 8.95 11.2 ...
##  $ V3 : num  55.7 58.2 56.9 53.7 56.5 50.9 57.8 45.7 48.2 56.3 ...
##  $ V4 : num  4.1 1.6 2.7 5.6 5.7 5.1 4.6 5.4 4.3 6.3 ...
##  $ V5 : num  9 3.8 8.1 18.9 34.5 21.9 16.7 60.5 24.4 29.6 ...
##  $ V6 : num  39.6 51.7 74 122.8 88.9 ...
##  $ V7 : int  279 80 107 147 180 150 186 640 182 85 ...
##  $ V8 : int  2 2 2 2 2 2 2 1 2 2 ...
##  $ V9 : int  4 2 3 4 1 2 3 2 3 1 ...
##  $ V10: int  207 51 82 53 134 147 151 399 130 59 ...
##  $ V11: int  241 52 54 148 151 106 129 360 118 66 ...
##  $ V12: num  60 40 20 40 40 40 40 60 40 40 ...
colnames(scenic) <- c("IDnumber","LengthStay","Age","InfectRisk","CultRatio","XrayRatio","NBeds","MedSchool","Region","DailyCensus","NNurses","Facilities")
names(scenic)
##  [1] "IDnumber"    "LengthStay"  "Age"         "InfectRisk"  "CultRatio"  
##  [6] "XrayRatio"   "NBeds"       "MedSchool"   "Region"      "DailyCensus"
## [11] "NNurses"     "Facilities"
#Summary para analisar quais variáveis terão que ser alteradas (de numéricas para categóricas, por exemplo)
summary(scenic)
##     IDnumber     LengthStay          Age          InfectRisk      CultRatio    
##  Min.   :  1   Min.   : 6.700   Min.   :38.80   Min.   :1.300   Min.   : 1.60  
##  1st Qu.: 29   1st Qu.: 8.340   1st Qu.:50.90   1st Qu.:3.700   1st Qu.: 8.40  
##  Median : 57   Median : 9.420   Median :53.20   Median :4.400   Median :14.10  
##  Mean   : 57   Mean   : 9.648   Mean   :53.23   Mean   :4.355   Mean   :15.79  
##  3rd Qu.: 85   3rd Qu.:10.470   3rd Qu.:56.20   3rd Qu.:5.200   3rd Qu.:20.30  
##  Max.   :113   Max.   :19.560   Max.   :65.90   Max.   :7.800   Max.   :60.50  
##    XrayRatio          NBeds         MedSchool        Region     
##  Min.   : 39.60   Min.   : 29.0   Min.   :1.00   Min.   :1.000  
##  1st Qu.: 69.50   1st Qu.:106.0   1st Qu.:2.00   1st Qu.:2.000  
##  Median : 82.30   Median :186.0   Median :2.00   Median :2.000  
##  Mean   : 81.63   Mean   :252.2   Mean   :1.85   Mean   :2.363  
##  3rd Qu.: 94.10   3rd Qu.:312.0   3rd Qu.:2.00   3rd Qu.:3.000  
##  Max.   :133.50   Max.   :835.0   Max.   :2.00   Max.   :4.000  
##   DailyCensus       NNurses        Facilities   
##  Min.   : 20.0   Min.   : 14.0   Min.   : 5.70  
##  1st Qu.: 68.0   1st Qu.: 66.0   1st Qu.:31.40  
##  Median :143.0   Median :132.0   Median :42.90  
##  Mean   :191.4   Mean   :173.2   Mean   :43.16  
##  3rd Qu.:252.0   3rd Qu.:218.0   3rd Qu.:54.30  
##  Max.   :791.0   Max.   :656.0   Max.   :80.00
#codificando corretamente as variáveis necessárias e gerando um novo resumo
scenic$IDnumber <- as.factor(scenic$IDnumber)
scenic$MedSchool <- as.factor(scenic$MedSchool)
scenic$Region <- as.factor(scenic$Region)

summary(scenic)
##     IDnumber     LengthStay          Age          InfectRisk      CultRatio    
##  1      :  1   Min.   : 6.700   Min.   :38.80   Min.   :1.300   Min.   : 1.60  
##  2      :  1   1st Qu.: 8.340   1st Qu.:50.90   1st Qu.:3.700   1st Qu.: 8.40  
##  3      :  1   Median : 9.420   Median :53.20   Median :4.400   Median :14.10  
##  4      :  1   Mean   : 9.648   Mean   :53.23   Mean   :4.355   Mean   :15.79  
##  5      :  1   3rd Qu.:10.470   3rd Qu.:56.20   3rd Qu.:5.200   3rd Qu.:20.30  
##  6      :  1   Max.   :19.560   Max.   :65.90   Max.   :7.800   Max.   :60.50  
##  (Other):107                                                                   
##    XrayRatio          NBeds       MedSchool Region  DailyCensus   
##  Min.   : 39.60   Min.   : 29.0   1:17      1:28   Min.   : 20.0  
##  1st Qu.: 69.50   1st Qu.:106.0   2:96      2:32   1st Qu.: 68.0  
##  Median : 82.30   Median :186.0             3:37   Median :143.0  
##  Mean   : 81.63   Mean   :252.2             4:16   Mean   :191.4  
##  3rd Qu.: 94.10   3rd Qu.:312.0                    3rd Qu.:252.0  
##  Max.   :133.50   Max.   :835.0                    Max.   :791.0  
##                                                                   
##     NNurses        Facilities   
##  Min.   : 14.0   Min.   : 5.70  
##  1st Qu.: 66.0   1st Qu.:31.40  
##  Median :132.0   Median :42.90  
##  Mean   :173.2   Mean   :43.16  
##  3rd Qu.:218.0   3rd Qu.:54.30  
##  Max.   :656.0   Max.   :80.00  
## 
#para utilizar os nomes das variáveis diretamente
attach(scenic)

#boxplot e histograma com linha de densidade estimada da variável de resposta
boxplot(LengthStay, horizontal = TRUE)

hist(LengthStay, freq = FALSE, breaks = "FD", main = "")
lines(density(LengthStay), col="blue", lwd=3)

#matriz de gráficos de dispersão
plot(scenic[, -c(1,8,9)])

#carregar biblioteca para fazer matriz de correlação
library(corrplot)
## corrplot 0.84 loaded
res <- cor(scenic[, -c(1,8,9)])
round(res,2)
##             LengthStay   Age InfectRisk CultRatio XrayRatio NBeds DailyCensus
## LengthStay        1.00  0.19       0.53      0.33      0.38  0.41        0.47
## Age               0.19  1.00       0.00     -0.23     -0.02 -0.06       -0.05
## InfectRisk        0.53  0.00       1.00      0.56      0.45  0.36        0.38
## CultRatio         0.33 -0.23       0.56      1.00      0.42  0.14        0.14
## XrayRatio         0.38 -0.02       0.45      0.42      1.00  0.05        0.06
## NBeds             0.41 -0.06       0.36      0.14      0.05  1.00        0.98
## DailyCensus       0.47 -0.05       0.38      0.14      0.06  0.98        1.00
## NNurses           0.34 -0.08       0.39      0.20      0.08  0.92        0.91
## Facilities        0.36 -0.04       0.41      0.19      0.11  0.79        0.78
##             NNurses Facilities
## LengthStay     0.34       0.36
## Age           -0.08      -0.04
## InfectRisk     0.39       0.41
## CultRatio      0.20       0.19
## XrayRatio      0.08       0.11
## NBeds          0.92       0.79
## DailyCensus    0.91       0.78
## NNurses        1.00       0.78
## Facilities     0.78       1.00
corrplot(res, method = "circle", order = "hclust", sig.level = 0.01, insig = "blank")

#A matriz obtida sugere correlação entre algumas variáveis, cujos diagramas de dispersão serão exibidos a seguir.
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
scatterplot(LengthStay ~ DailyCensus|MedSchool, data = scenic)

scatterplot(LengthStay ~ InfectRisk|MedSchool, data = scenic)

scatterplot(LengthStay ~ Facilities|MedSchool, data = scenic)

#analisando como a variável de resposta (LengthStay) depende das outras variáveis

#comparando a variável de resposta com outras variáveis

boxplot(LengthStay ~ InfectRisk, horizontal = TRUE)

boxplot(LengthStay ~ Facilities, horizontal = TRUE)

#Diagrama de dispersão condicional para cada região
scatterplot(LengthStay ~ CultRatio|Region, data = scenic)

scatterplot(LengthStay ~ XrayRatio|Region, data = scenic)

scatterplot(LengthStay ~ Facilities|Region, data = scenic)

scatterplot(LengthStay ~ InfectRisk|Region, data = scenic)

#Construção de nova base de dados

data_rls <- data.frame(LengthStay = LengthStay, InfectRisk = InfectRisk,Facilities = Facilities,XrayRatio = XrayRatio )
summary(data_rls)
##    LengthStay       InfectRisk      Facilities      XrayRatio     
##  Min.   : 6.700   Min.   :1.300   Min.   : 5.70   Min.   : 39.60  
##  1st Qu.: 8.340   1st Qu.:3.700   1st Qu.:31.40   1st Qu.: 69.50  
##  Median : 9.420   Median :4.400   Median :42.90   Median : 82.30  
##  Mean   : 9.648   Mean   :4.355   Mean   :43.16   Mean   : 81.63  
##  3rd Qu.:10.470   3rd Qu.:5.200   3rd Qu.:54.30   3rd Qu.: 94.10  
##  Max.   :19.560   Max.   :7.800   Max.   :80.00   Max.   :133.50
# Constroi modelo de regressao linear simples (rls)

#Para InfectRisk

ir_rls <- lm(LengthStay ~ InfectRisk, data = data_rls)
ir_rls
## 
## Call:
## lm(formula = LengthStay ~ InfectRisk, data = data_rls)
## 
## Coefficients:
## (Intercept)   InfectRisk  
##      6.3368       0.7604
summary(ir_rls)
## 
## Call:
## lm(formula = LengthStay ~ InfectRisk, data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0587 -0.7776 -0.1487  0.7159  8.2805 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.3368     0.5213  12.156  < 2e-16 ***
## InfectRisk    0.7604     0.1144   6.645 1.18e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.624 on 111 degrees of freedom
## Multiple R-squared:  0.2846, Adjusted R-squared:  0.2781 
## F-statistic: 44.15 on 1 and 111 DF,  p-value: 1.177e-09
anova(ir_rls)
# Adiciona curva de regressao estimada no grafico de dispersao
plot(LengthStay ~InfectRisk)
abline(ir_rls$coef, col=2, lwd=3)

#Para Facilities
fa_rls <- lm(LengthStay ~ Facilities, data = data_rls)
fa_rls
## 
## Call:
## lm(formula = LengthStay ~ Facilities, data = data_rls)
## 
## Coefficients:
## (Intercept)   Facilities  
##     7.71877      0.04471
summary(fa_rls)
## 
## Call:
## lm(formula = LengthStay ~ Facilities, data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2712 -1.0716 -0.2816  0.7584  9.5433 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.71877    0.51020  15.129  < 2e-16 ***
## Facilities   0.04471    0.01116   4.008 0.000111 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.795 on 111 degrees of freedom
## Multiple R-squared:  0.1264, Adjusted R-squared:  0.1185 
## F-statistic: 16.06 on 1 and 111 DF,  p-value: 0.0001113
anova(fa_rls)
# Adiciona curva de regressao estimada no grafico de dispersao
plot(LengthStay ~Facilities)
abline(fa_rls$coef, col=2, lwd=3)

#Para XrayRatio
xr_rls <- lm(LengthStay ~ XrayRatio, data = data_rls)
xr_rls
## 
## Call:
## lm(formula = LengthStay ~ XrayRatio, data = data_rls)
## 
## Coefficients:
## (Intercept)    XrayRatio  
##     6.56637      0.03776
summary(xr_rls)
## 
## Call:
## lm(formula = LengthStay ~ XrayRatio, data = data_rls)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9226 -1.0810 -0.2708  0.8200  8.7008 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.566373   0.726094   9.043 5.67e-15 ***
## XrayRatio   0.037756   0.008657   4.361 2.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.774 on 111 degrees of freedom
## Multiple R-squared:  0.1463, Adjusted R-squared:  0.1386 
## F-statistic: 19.02 on 1 and 111 DF,  p-value: 2.906e-05
anova(xr_rls)
# Adiciona curva de regressao estimada no grafico de dispersao
plot(LengthStay ~XrayRatio)
abline(xr_rls$coef, col=2, lwd=3)

#MSE
#os valores de MSE encontrados:
anova(ir_rls)[1,"Mean Sq"]
## [1] 116.4459
anova(fa_rls)[1,"Mean Sq"]
## [1] 51.72714
anova(xr_rls)[1,"Mean Sq"]
## [1] 59.86438

Com isso, concluimos que a variavel explicativa que produz menor variabilidade entorno da reta de regressão é a Facilities.

Analisando os valores de R^2, nenhuma variavél possui contribuição satisfatória para explicar o LengthStay, entretanto a com melhor coeficiente de determinação é o InfectRisk.

data_rls <- data.frame(LengthStay = LengthStay, InfectRisk = InfectRisk, Region = Region)
summary(data_rls)
##    LengthStay       InfectRisk    Region
##  Min.   : 6.700   Min.   :1.300   1:28  
##  1st Qu.: 8.340   1st Qu.:3.700   2:32  
##  Median : 9.420   Median :4.400   3:37  
##  Mean   : 9.648   Mean   :4.355   4:16  
##  3rd Qu.:10.470   3rd Qu.:5.200         
##  Max.   :19.560   Max.   :7.800
attach(data_rls)
## The following objects are masked from scenic:
## 
##     InfectRisk, LengthStay, Region
datareg1 <- subset(data_rls, Region == 1)
datareg2 <- subset(data_rls, Region == 2)
datareg3 <- subset(data_rls, Region == 3)
datareg4 <- subset(data_rls, Region == 4)

scenic_rls1 <-lm(datareg1$LengthStay ~ datareg1$InfectRisk)
scenic_rls2 <-lm(datareg2$LengthStay ~ datareg2$InfectRisk)
scenic_rls3 <-lm(datareg3$LengthStay ~ datareg3$InfectRisk)
scenic_rls4 <-lm(datareg4$LengthStay ~ datareg4$InfectRisk)
plot(LengthStay ~ InfectRisk,  xlim = c(0,8), ylim = c(6,20) )
legend("topleft", legend = c("Região 1","Região 2", "Região 3", "Região 4"), col = c("red", "blue", "green", "yellow"), lty = 1, cex = 0.8)
abline(scenic_rls1, col="red", lwd=3)
abline(scenic_rls2, col="blue", lwd=3)
abline(scenic_rls3, col="green", lwd=3)
abline(scenic_rls4, col="yellow", lwd=3)

Vamos separar os gráficos em 4, para poder fazer uma análise mais precisa.

plot(datareg1$LengthStay ~ datareg1$InfectRisk, xlim = c(0,8), ylim = c(6,20), xlab = "InfectRisk", ylab = "LengthStay")
legend("topleft", legend = c("Região 1"))
abline(scenic_rls1, col="red", lwd=3)

plot(datareg2$LengthStay ~ datareg2$InfectRisk, xlim = c(0,8), ylim = c(6,20), xlab = "InfectRisk", ylab = "LengthStay" )
legend("topleft", legend = c("Região 2"))
abline(scenic_rls2, col="blue", lwd=3)

plot(datareg3$LengthStay ~ datareg3$InfectRisk, xlim = c(0,8), ylim = c(6,20), xlab = "InfectRisk", ylab = "LengthStay" )
legend("topleft", legend = c("Região 3"))
abline(scenic_rls3, col="green", lwd=3)

plot(datareg4$LengthStay ~ datareg4$InfectRisk, xlim = c(0,8), ylim = c(6,20), xlab = "InfectRisk", ylab = "LengthStay" )
legend("topleft", legend = c("Região 4"))
abline(scenic_rls4, col="yellow", lwd=3)

De acordo com o primeiro gráfico, as funções de regressão foram semelhantes nas regiões 2 e 3 e, nas regiões 1 e 4, foram diferentes entre si e das demais. Analisando as retas, sugere-se que o InfectRisk não tem impacto no LengthStay para a região 4, enquanto que, para a região 1, a relação é direta, visto que o aumento do InfectRisk faz crescer o valor de LengthStay; já para as regiões 2 e 3, nota-se uma possível relação, mas mais sutil que para a Região 1. Contudo, é necessário se atentar a outros fatores relevantes para garantir a qualidade da regressão e a veracidade da correlação das variáveis.Tal análise será feita a seguir.

anova(scenic_rls1)[1,"Mean Sq"]
## [1] 79.24452
anova(scenic_rls2)[1,"Mean Sq"]
## [1] 12.97975
anova(scenic_rls3)[1,"Mean Sq"]
## [1] 21.12862
anova(scenic_rls4)[1,"Mean Sq"]
## [1] 0.0034406

A variabilidade em torno das retas são diferentes para cada região. Nota-se que, para a Regiâo 4, a variabilidade é baixa, o que sustenta a hipótese de que o InfectRisk não tem impacto no LengthStay para essa Região. Devido à alta variabilidade em torno da reta da Região 1, não se pode afirma que a alta inclinação implica em alta influência do crescimento do LengthStay em função do crescimento do InfectRisk. Já para as Regiões 2 e 3, o valor do MSE foi baixo em relação ao da Região 1, o que indica melhor qualidade na regressão, entretanto, os valores, ainda assim, foram bem mais altos do que o MSE da Região 4.

Faremos “summary” para cada subconjunto de dados respectivo a cada Região para conferir os valores dois coeficientes em relação aos intervalos de confiança.

summary(scenic_rls1)
## 
## Call:
## lm(formula = datareg1$LengthStay ~ datareg1$InfectRisk)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1887 -1.0713 -0.3449  0.6822  6.2617 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           4.5379     1.5852   2.863 0.008193 ** 
## datareg1$InfectRisk   1.3477     0.3159   4.267 0.000233 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.086 on 26 degrees of freedom
## Multiple R-squared:  0.4118, Adjusted R-squared:  0.3892 
## F-statistic:  18.2 on 1 and 26 DF,  p-value: 0.0002326
summary(scenic_rls2)
## 
## Call:
## lm(formula = datareg2$LengthStay ~ datareg2$InfectRisk)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1998 -0.6047  0.1244  0.7435  1.8283 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           7.5605     0.6267  12.063 4.89e-13 ***
## datareg2$InfectRisk   0.4832     0.1366   3.536  0.00134 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.019 on 30 degrees of freedom
## Multiple R-squared:  0.2942, Adjusted R-squared:  0.2707 
## F-statistic: 12.51 on 1 and 30 DF,  p-value: 0.001341
summary(scenic_rls3)
## 
## Call:
## lm(formula = datareg3$LengthStay ~ datareg3$InfectRisk)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9748 -0.4921 -0.2071  0.2900  2.4954 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           7.1293     0.4632  15.393  < 2e-16 ***
## datareg3$InfectRisk   0.5251     0.1107   4.742 3.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9694 on 35 degrees of freedom
## Multiple R-squared:  0.3911, Adjusted R-squared:  0.3737 
## F-statistic: 22.48 on 1 and 35 DF,  p-value: 3.494e-05
summary(scenic_rls4)
## 
## Call:
## lm(formula = datareg4$LengthStay ~ datareg4$InfectRisk)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4158 -0.6716 -0.3077  0.6918  2.0425 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          8.03805    1.36481   5.889 3.94e-05 ***
## datareg4$InfectRisk  0.01728    0.30583   0.056    0.956    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.038 on 14 degrees of freedom
## Multiple R-squared:  0.0002279,  Adjusted R-squared:  -0.07118 
## F-statistic: 0.003192 on 1 and 14 DF,  p-value: 0.9557
confint.lm(scenic_rls1, level=0.95)
##                         2.5 %   97.5 %
## (Intercept)         1.2794501 7.796392
## datareg1$InfectRisk 0.6984451 1.997047
confint.lm(scenic_rls2, level=0.95)
##                         2.5 %    97.5 %
## (Intercept)         6.2805261 8.8404863
## datareg2$InfectRisk 0.2041384 0.7622031
confint.lm(scenic_rls3, level=0.95)
##                         2.5 %   97.5 %
## (Intercept)         6.1890595 8.069615
## datareg3$InfectRisk 0.3002664 0.749899
confint.lm(scenic_rls4, level=0.95)
##                          2.5 %     97.5 %
## (Intercept)          5.1108117 10.9652846
## datareg4$InfectRisk -0.6386564  0.6732136

As retas de regressão para cada região não parecem ter mesma inclinação, devido a IC’s discrepantes. Podemos concluir que, em cada Região, o impacto do valor de InfectRisk é diferente sobre o valor de LengthStay.

#criação de novos lm para que a criação dos intervalos de confiança funcione
scenic_reg1 <- lm(LengthStay ~ InfectRisk, data = datareg1)
scenic_reg2 <- lm(LengthStay ~ InfectRisk, data = datareg2)
scenic_reg3 <- lm(LengthStay ~ InfectRisk, data = datareg3)
scenic_reg4 <- lm(LengthStay ~ InfectRisk, data = datareg4)

Xc1 <- data.frame(InfectRisk = 5)
predict.lm(scenic_reg1, newdata=Xc1, interval="confidence", level = 0.95)
##        fit      lwr      upr
## 1 11.27665 10.46114 12.09216
Xc2 <- data.frame(InfectRisk = 5)
predict.lm(scenic_reg2, newdata=Xc2, interval="confidence", level = 0.95)
##       fit      lwr     upr
## 1 9.97636 9.571522 10.3812
Xc3 <- data.frame(InfectRisk = 5)
predict.lm(scenic_reg3, newdata=Xc3, interval="confidence", level = 0.95)
##        fit     lwr      upr
## 1 9.754751 9.35118 10.15832
Xc4 <- data.frame(InfectRisk = 5)
predict.lm(scenic_reg4, newdata=Xc4, interval="confidence", level = 0.95)
##        fit      lwr      upr
## 1 8.124441 7.435514 8.813368

Concluímos que com 95% de confiança que, com ‘InfectRisk = 5’, na:

Região 1, o valor de LengthStay está no intervalo de 10.46114 a 12.09216 dias;

Região 2, o valor de LengthStay está no intervalo de 9.571522 a 10.3812 dias;

Região 3, o valor de LengthStay está no intervalo de 9.35118 a 10.15832 dias;

Região 4, o valor de LengthStay está no intervalo de 7.435514 a 8.813368 dias.

predict.lm(scenic_reg1, newdata = Xc1, interval = "prediction", level = 0.95)
##        fit      lwr      upr
## 1 11.27665 6.911133 15.64217
predict.lm(scenic_reg2, newdata = Xc2, interval = "prediction", level = 0.95)
##       fit      lwr      upr
## 1 9.97636 7.856747 12.09597
predict.lm(scenic_reg3, newdata = Xc3, interval = "prediction", level = 0.95)
##        fit      lwr      upr
## 1 9.754751 7.745751 11.76375
predict.lm(scenic_reg4, newdata = Xc4, interval = "prediction", level = 0.95)
##        fit      lwr      upr
## 1 8.124441 5.793563 10.45532

Concluímos que com 95% de confiança que, com ‘InfectRisk = 5’, na:

Região 1, o LengthStay esperado para um novo hospital está eno intervalo de 6.911133 a 15.64217 dias;

Região 2, o LengthStay esperado para um novo hospital está eno intervalo de 7.856747 a 12.095977 dias;

Região 3, o LengthStay esperado para um novo hospital está eno intervalo de 7.745751 a 11.76375 dias;

Região 4, o LengthStay esperado para um novo hospital está eno intervalo de 5.793563 a 10.45532 dias.

Diagnóstico

Diagnóstico através da análise dos resíduos para InfectRisk

par(mfcol = c(2,2), mar=c(4,4,3,1))

#Gráficos dos resíduos:

#para dados originais: resid x lengthstay_orig; resid x infectrisk
resid_rls_orig <-  ir_rls$res
lengthstay_orig <- ir_rls$fit

plot(resid_rls_orig ~ lengthstay_orig, main = "dados originais", ylab="resíduos")
abline(h = 0, col = "red", lty = "dashed")
plot(resid_rls_orig ~ InfectRisk, ylab ="resíduos")
abline(h = 0, col = "red", lty = "dashed")

Diagnóstico através da análise dos resíduos para Facilities

par(mfcol = c(2,2), mar=c(4,4,3,1))

#Gráficos dos resíduos:

#para dados originais: resid x lengthstay_orig; resid x facilities
resid_rls_orig <-  fa_rls$res
lengthstay_orig <- fa_rls$fit

plot(resid_rls_orig ~ lengthstay_orig, main = "dados originais", ylab="resíduos")
abline(h = 0, col = "red", lty = "dashed")
plot(resid_rls_orig ~ Facilities, ylab ="resíduos")
abline(h = 0, col = "red", lty = "dashed")

Diagnóstico através da análise dos resíduos para XrayRatio

par(mfcol = c(2,2), mar=c(4,4,3,1))

#Gráficos dos resíduos:

#para dados originais: resid x lengthstay_orig; resid x xrayratio
resid_rls_orig <-  xr_rls$res
lengthstay_orig <- xr_rls$fit

plot(resid_rls_orig ~ lengthstay_orig, main = "dados originais", ylab="resíduos")
abline(h = 0, col = "red", lty = "dashed")
plot(resid_rls_orig ~ XrayRatio, ylab ="resíduos")
abline(h = 0, col = "red", lty = "dashed")

Note que os gráficos de resíduos para os dados originais versus variável de resposta ajustada apresentam padrão linear nas três variáveis explicativas(InfectRisk, Facilities e XrayRatio).Tal fato sugere que não é necessário fazer uma transformação dos dados e consideraremos, portanto, o modelo original.

Vamos fazer os gráficos do módulo dos resíduos para cada uma das três variáveis explicativas para permitir analisar efeitos com mais facilidade.

InfectRisk:

par(mfcol = c(2,2), mar=c(4,4,3,1))

#Gráficos dos resíduos:

#para dados originais: resid x lengthstay_orig; resid x infectrisk
resid_rls_orig <-  ir_rls$res
lengthstay_orig <- ir_rls$fit

plot(abs(resid_rls_orig) ~ lengthstay_orig, main = "dados originais", ylab="|resíduos|")
abline(h = 0, col = "red", lty = "dashed")
plot(abs(resid_rls_orig) ~ InfectRisk, ylab ="|resíduos|")
abline(h = 0, col = "red", lty = "dashed")

Facilities:

par(mfcol = c(2,2), mar=c(4,4,3,1))

#Gráficos dos resíduos:

#para dados originais: resid x lengthstay_orig; resid x facilities
resid_rls_orig <-  fa_rls$res
lengthstay_orig <- fa_rls$fit

plot(abs(resid_rls_orig) ~ lengthstay_orig, main = "dados originais", ylab="|resíduos|")
abline(h = 0, col = "red", lty = "dashed")
plot(abs(resid_rls_orig) ~ Facilities, ylab ="|resíduos|")
abline(h = 0, col = "red", lty = "dashed")

XrayRatio:

par(mfcol = c(2,2), mar=c(4,4,3,1))

#Gráficos dos resíduos:

#para dados originais: resid x lengthstay_orig; resid x xrayratio
resid_rls_orig <-  xr_rls$res
lengthstay_orig <- xr_rls$fit

plot(abs(resid_rls_orig) ~ lengthstay_orig, main = "dados originais", ylab="|resíduos|")
abline(h = 0, col = "red", lty = "dashed")
plot(abs(resid_rls_orig) ~ XrayRatio, ylab ="|resíduos|")
abline(h = 0, col = "red", lty = "dashed")

Os gráficos para as três variáveis sugerem que a variância está aproximadamente constante.Para verificar formalmente, vamos fazer o Teste de Breusch-Pagan com hipótese nula de que a variância é constante; a hipótese alternativa é de que a variância não se mantém constante para todas as observações.

Teste de Homoscedasticidade de Breusch-Pagan

Ho: sigma^2 = cte

Ha: sigma^2 != cte

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# InfectRisk

bptest(ir_rls)
## 
##  studentized Breusch-Pagan test
## 
## data:  ir_rls
## BP = 4.8996, df = 1, p-value = 0.02686
#Facilities

bptest(fa_rls)
## 
##  studentized Breusch-Pagan test
## 
## data:  fa_rls
## BP = 2.1424, df = 1, p-value = 0.1433
#XrayRatio

bptest(xr_rls)
## 
##  studentized Breusch-Pagan test
## 
## data:  xr_rls
## BP = 3.883, df = 1, p-value = 0.04878

Os p-valores implicam que, para as variáveis InfectRisk e XrayRatio, a hipótese nula é descartada e, para a variável Facilities, a hipótese nula é aceita, sob um fator de risco de 5%. Logo, conclui-se que para InfectRisk e XrayRatio, sugere-se que a variância não é constante. Já para Facilities, é provável que seja constante.

#Outliers Para identificar possíveis outliers, vamos analisar os resíduos padronizados (os resíduos semi-studentizados produzem praticamente os mesmos gráficos):

InfectRisk:

par(mfrow=c(1,2))

# Gráficos para os resíduos padronizados
resid_rls_pad_ir <- rstandard(ir_rls)

plot(resid_rls_pad_ir ~ lengthstay_orig, ylim=c(-5,5),
     ylab="resíduos padronizados")
abline(h = c(-4, -3, 0, 3, 4) , col = "red", lty = "dashed")

Facilities:

par(mfrow=c(1,2))

# Gráficos para os resíduos padronizados
resid_rls_pad_f <- rstandard(fa_rls)

plot(resid_rls_pad_f ~ lengthstay_orig, ylim=c(-5,5),
     ylab="resíduos padronizados")
abline(h = c(-4, -3, 0, 3, 4) , col = "red", lty = "dashed")

XrayRatio:

par(mfrow=c(1,2))

# Gráficos para os resíduos padronizados
resid_rls_pad_xr <- rstandard(xr_rls)

plot(resid_rls_pad_xr ~ lengthstay_orig, ylim=c(-5,5),
     ylab="resíduos padronizados")
abline(h = c(-4, -3, 0, 3, 4) , col = "red", lty = "dashed")

Há algumas observações extremas (possíveis outliers), mas elas não devem ser eliminadas neste ponto da análise. Antes disso, é necessário avaliar o efeito da inclusão de outras variáveis explicativas no modelo.

#Independência

Assumimos que os erros não são correlacionados. Não é possível identificar padrões aparentes sugestivos de existência de correlação entre os erros. O procedimento formal para avaliar esta hipótese é o Teste de Durbin-Watson (a hipótese nula é de que a correlação é nula e a hipótese alternativa automática é de que a correlação é positiva, mas nós realizaremos o teste bi-caudal):

# Teste de Durbin-Watson para correlação nula dos erros
# Ho: corr  = 0
# Ha: corr != 0

library(lmtest)

#InfectRisk
dwtest(ir_rls, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  ir_rls
## DW = 2.0964, p-value = 0.6129
## alternative hypothesis: true autocorrelation is not 0
#Facilities
dwtest(fa_rls, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  fa_rls
## DW = 2.0521, p-value = 0.7818
## alternative hypothesis: true autocorrelation is not 0
#XrayRatio
dwtest(xr_rls, alternative = "two.sided")
## 
##  Durbin-Watson test
## 
## data:  xr_rls
## DW = 1.9201, p-value = 0.663
## alternative hypothesis: true autocorrelation is not 0

Portanto, devido aos valores-p encontrados, aceita-se a hipótese nula para todas as três variáveis, ou seja, sugere-se que os erros não são correlacionados.

#Normalidade

Podemos avaliar a normalidade através de métodos gráficos e testes de hipóteses. O procedimento gráfico se baseia na análise de histogramas, boxplots e gráficos de quantis para os resíduos (padronizados).

InfectRisk:

par(mfrow=c(1,2))

# Histograma
hist(resid_rls_pad_ir)

# Gráfico de quantis
qqnorm(resid_rls_pad_ir)
qqline(resid_rls_pad_ir, lty="dashed")

Facilities:

par(mfrow=c(1,2))

# Histograma
hist(resid_rls_pad_f)

# Gráfico de quantis
qqnorm(resid_rls_pad_f)
qqline(resid_rls_pad_f, lty="dashed")

XrayRatio:

par(mfrow=c(1,2))

# Histograma
hist(resid_rls_pad_xr)

# Gráfico de quantis
qqnorm(resid_rls_pad_xr)
qqline(resid_rls_pad_xr, lty="dashed")

#Fazer comentário do histograma e do qqplot

Visualmente, o histograma não aparenta ter um caráter normal de distribuição. A hipótese de normalidade pode ser formalmente testada através do teste de Shapiro-Wilk. A hipótese nula é de que os resíduos são normalmente distribuídos, versus a hipótese alternativa de que os resíduos não seguem a distribuição normal:

# Teste de Normalidade de Shapiro-Wilk
# Ho: normal
# Ha: não-normal

shapiro.test(resid_rls_pad_ir)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_rls_pad_ir
## W = 0.86898, p-value = 1.459e-08
shapiro.test(resid_rls_pad_f)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_rls_pad_f
## W = 0.8615, p-value = 7.129e-09
shapiro.test(resid_rls_pad_xr)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_rls_pad_xr
## W = 0.87457, p-value = 2.534e-08

O valor-p do teste nos leva a rejeitar a hipótese nula e concluir que os erros não são normalmente distribuídos. Como sabemos, a violação da hipótese de normalidade não é das mais graves, contudo, sugere-se a necessidade de outras medidas corretivas.

Faremos agora a tranformação Box-Cox nessa base de dados para testar a linearidade do modelo.

Transformação Box-Cox

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
 #InfectRIsk
bcir <- boxcox(lm(LengthStay ~ InfectRisk), lambda = seq(-2, 2, by = 0.1), plotit = TRUE)

(lambdair <- bcir$x[which.max(bcir$y)])
## [1] -1.111111
#Facilities
bcfa <- boxcox(lm(LengthStay ~ Facilities), lambda = seq(-2, 2, by = 0.1), plotit = TRUE)

(lambdafa <- bcfa$x[which.max(bcfa$y)])
## [1] -1.353535
#XrayRatio  
bcxr <- boxcox(lm(LengthStay ~ XrayRatio), lambda = seq(-2, 2, by = 0.1), plotit = TRUE)

(lambdaxr <- bcxr$x[which.max(bcxr$y)])
## [1] -1.232323

InfectRisk: O valor de \(\lambda\) que maximiza a função de verossimilhança (ou melhor, o logaritmo dessa função) é próximo de -1.1, o que indica que a transformação \(\frac{1}{Y^{1.1}}\) deve ser aplicada. Sendo assim, replicaremos o procedimento realizado para as variáveis transformadas:

#Transforma da variável de resposta

LengthStay_bcir <- 1/(LengthStay^(1.1))

# Constrói nova base de dados
data_rls_bcir <- data.frame(InfectRisk = InfectRisk, LengthStay_bcir = LengthStay_bcir )

# Ajusta modelo de regressão linear simples aos dados transformados
scenic_rls_bcir <- lm(LengthStay_bcir ~ InfectRisk, data = data_rls_bcir)

# Gera resumo do modelo de regressão ajustado
summary(scenic_rls_bcir)
## 
## Call:
## lm(formula = LengthStay_bcir ~ InfectRisk, data = data_rls_bcir)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.034087 -0.007782 -0.000250  0.006908  0.038646 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.113375   0.004241  26.734  < 2e-16 ***
## InfectRisk  -0.006360   0.000931  -6.831 4.74e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01321 on 111 degrees of freedom
## Multiple R-squared:  0.296,  Adjusted R-squared:  0.2896 
## F-statistic: 46.66 on 1 and 111 DF,  p-value: 4.736e-10

Modelo ajustado:

\(\frac{1}{LengthStay^{1.1}} = -0.00636*InfectRisk + 0.113375\)

# Adiciona curva de regressão estimada ao gráfico de dispersão
plot(data_rls_bcir)
abline(scenic_rls_bcir$coef, col = 2, lwd = 3)

Precisamos realizar o diagnóstico do novo modelo:

par(mfrow = c(2,2), mar=c(4,4,3,1))
# Produz gráficos dos resíduos
plot(scenic_rls_bcir)

Facilities:

O valor de \(\lambda\) que maximiza a função de verossimilhança (ou melhor, o logaritmo dessa função) é próximo de -1.3, o que indica que a transformação \(\frac{1}{Y^{1.3}}\) deve ser aplicada. Sendo assim, replicaremos o procedimento realizado para as variáveis transformadas:

#Transforma da variável de resposta

LengthStay_bcfa <- 1/(LengthStay^(1.3))

# Constrói nova base de dados
data_rls_bcfa <- data.frame(Facilities = Facilities, LengthStay_bcfa = LengthStay_bcfa)

# Ajusta modelo de regressão linear simples aos dados transformados
scenic_rls_bcfa <- lm(LengthStay_bcfa ~ Facilities, data = data_rls_bcfa)

# Gera resumo do modelo de regressão ajustado
summary(scenic_rls_bcfa)
## 
## Call:
## lm(formula = LengthStay_bcfa ~ Facilities, data = data_rls_bcfa)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0315774 -0.0069620 -0.0000135  0.0072264  0.0278579 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.799e-02  3.106e-03  21.891  < 2e-16 ***
## Facilities  -3.007e-04  6.791e-05  -4.429 2.23e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01092 on 111 degrees of freedom
## Multiple R-squared:  0.1502, Adjusted R-squared:  0.1425 
## F-statistic: 19.61 on 1 and 111 DF,  p-value: 2.231e-05

Modelo ajustado:

\(\frac{1}{LengthStay^{1.3}} = -0.0003*Facilities + 0.068\)

# Adiciona curva de regressão estimada ao gráfico de dispersão
plot(data_rls_bcfa)
abline(scenic_rls_bcfa$coef, col = 2, lwd = 3)

Precisamos realizar o diagnóstico do novo modelo:

par(mfrow = c(2,2), mar=c(4,4,3,1))
# Produz gráficos dos resíduos
plot(scenic_rls_bcfa)

XrayRatio: O valor de \(\lambda\) que maximiza a função de verossimilhança (ou melhor, o logaritmo dessa função) é próximo de -1.2, o que indica que a transformação \(\frac{1}{Y^{1.2}}\) deve ser aplicada. Sendo assim, replicaremos o procedimento realizado para as variáveis transformadas:

#Transforma da variável de resposta

LengthStay_bcxr <- 1/(LengthStay^(1.2))

# Constrói nova base de dados
data_rls_bcxr <- data.frame(XrayRatio = XrayRatio, LengthStay_bcxr = LengthStay_bcxr)

# Ajusta modelo de regressão linear simples aos dados transformados
scenic_rls_bcxr <- lm(LengthStay_bcxr ~ XrayRatio, data = data_rls_bcxr)

# Gera resumo do modelo de regressão ajustado
summary(scenic_rls_bcxr)
## 
## Call:
## lm(formula = LengthStay_bcxr ~ XrayRatio, data = data_rls_bcxr)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.034671 -0.008227 -0.000078  0.007840  0.033166 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.015e-02  5.202e-03  17.331  < 2e-16 ***
## XrayRatio   -2.635e-04  6.202e-05  -4.249 4.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01271 on 111 degrees of freedom
## Multiple R-squared:  0.1399, Adjusted R-squared:  0.1321 
## F-statistic: 18.05 on 1 and 111 DF,  p-value: 4.489e-05

Modelo ajustado:

\(\frac{1}{LengthStay^{1.2}} = -0.00026*XrayRatio + 0.09\)

# Adiciona curva de regressão estimada ao gráfico de dispersão
plot(data_rls_bcxr)
abline(scenic_rls_bcxr$coef, col = 2, lwd = 3)

Precisamos realizar o diagnóstico do novo modelo:

par(mfrow = c(2,2), mar=c(4,4,3,1))
# Produz gráficos dos resíduos
plot(scenic_rls_bcxr)

Após fazer a transformação, podemos observar pelo gráfico dos resíduos que esta não apresenta nenhuma curvatura, indicando, assim, uma relação linear, melhorando o modelo.

##Criação da nova base de dados para remover outliers

newdata_rls <- data_rls[-c(47,112), ]
newir_rls <- lm(LengthStay ~ InfectRisk, data = newdata_rls)  
newlengthstay_orig <- newir_rls$fit  
par(mfrow=c(1,2))

# Gráficos para os resíduos padronizados
newresid_rls_pad <- rstandard(newir_rls)

plot(newresid_rls_pad ~ newlengthstay_orig, ylim=c(-5,5),
     ylab="resíduos padronizados")
abline(h = c(-4, -3, 0, 3, 4) , col = "red", lty = "dashed")

Construiremos intervalos de previsão de 95% de confiança para novas observações que apresentam valores da variável explicativa iguais a 6.5 e 5.9.

Xi1 <- data.frame(InfectRisk = 6.5)
predict.lm(newir_rls, newdata = Xi1, interval="confidence", level = 0.95)
##        fit      lwr      upr
## 1 10.81259 10.36399 11.26118
Xi2 <- data.frame(InfectRisk = 5.9)
predict.lm(newir_rls, newdata = Xi2, interval="confidence", level = 0.95)
##        fit      lwr      upr
## 1 10.44674 10.08424 10.80923

Concluimos que, com 95% de confiança, para InfectRisk = 6.5 o valor de LengthStay está entre 10.36399 e 11.26118 dias. Já para InfectRisk = 5.9, com 95% de confiança que LengthStay está entre 10.08424 e 10.80923 dias. Portanto, as observações removidas não estão no intervalo de confiança encontrado, pois eram de LengthStay = 19.56 e 17.94 dias. Logo, foi coerente a remoção de ambas observações.

Diagnóstico por Região

## Região 1

par(mfrow=c(1,1))

resid_rls_orig_reg1 <- scenic_rls1$res
lengthstay_orig_reg1 <- scenic_rls1$fit
plot(abs(resid_rls_orig_reg1) ~ lengthstay_orig_reg1, main = "Região 1", ylab="|resíduos|", xlim = c(7.5, 15), ylim = c(0, 7))
abline(h = 0, col = "red", lty="dashed")

## Região 2
par(mfrow=c(1,1))

resid_rls_orig_reg2 <- scenic_rls2$res
lengthstay_orig_reg2 <- scenic_rls2$fit
plot(abs(resid_rls_orig_reg2) ~ lengthstay_orig_reg2, main = "Região 2", ylab="|resíduos|", xlim = c(7.5, 15), ylim = c(0, 7))
abline(h = 0, col = "red", lty="dashed")

## Região 3
par(mfrow=c(1,1))

resid_rls_orig_reg3 <- scenic_rls3$res
lengthstay_orig_reg3 <- scenic_rls3$fit
plot(abs(resid_rls_orig_reg3) ~ lengthstay_orig_reg3, main = "Região 3", ylab="|resíduos|", xlim = c(7.5, 15), ylim = c(0, 7))
abline(h = 0, col = "red", lty="dashed")

## Região 4
par(mfrow=c(1,1))

resid_rls_orig_reg4 <- scenic_rls4$res
lengthstay_orig_reg4 <- scenic_rls4$fit
plot(abs(resid_rls_orig_reg4) ~ lengthstay_orig_reg4, main = "Região 4", ylab="|resíduos|", xlim = c(7.5, 15), ylim = c(0, 7))
abline(h = 0, col = "red", lty="dashed")

Teste de Homoscedasticidade de Breusch-Pagan para cada Região:

## Ho: sigma^2 = cte
## Ha: sigma^2 != cte
## Região 1:
bptest(scenic_rls1)
## 
##  studentized Breusch-Pagan test
## 
## data:  scenic_rls1
## BP = 4.6187, df = 1, p-value = 0.03163
## Região 2:
bptest(scenic_rls2)
## 
##  studentized Breusch-Pagan test
## 
## data:  scenic_rls2
## BP = 0.21, df = 1, p-value = 0.6468
## Região 3:
bptest(scenic_rls3)
## 
##  studentized Breusch-Pagan test
## 
## data:  scenic_rls3
## BP = 0.0083117, df = 1, p-value = 0.9274
## Região 4:
bptest(scenic_rls4)
## 
##  studentized Breusch-Pagan test
## 
## data:  scenic_rls4
## BP = 2.1159, df = 1, p-value = 0.1458

Os p-valores implicam que, para a Região 1, a hipótese nula é descartada e, para as Regiões 2, 3 e 4, a hipótese nula é aceita, sob um fator de risco de 5%. Logo, conclui-se que para a Região 1, sugere-se que a variância não é constante. Já para as Regiões 2, 3 e 4, é provável que seja constante. Visualmente, pelos gráficos do módulo dos resíduos pelos valores de InfectRisk, percebemos que apenas a Região 1 sugere heterocedasticidade, e o contrário para as demais, o que corrobora com o resultado do teste realizado.

Análise

  1. Assuma que um modelo de regressão linear simples é adequado para modelar a relação da variável de resposta LengthStay a cada uma das variáveis explicativas InfectRisk, Facilities e XrayRatio.
  • Construa um modelo de regressão para cada um desses pares de variáveis;
  • Construa gráficos de dispersão (separados) com as retas de regressão ajustadas para cada caso;
  • Calcule o MSE para cada modelo. Que variável explicativa produz menor variabilidade em torno da reta de regressão ajustada?
  • Utilizando R2 como critério, qual das variáveis explicativas contribui para a maior redução na variabilidade da resposta LengthStay?
  1. Para cada região geográfica, construa um modelo de regressão para a variável de resposta LengthStay em função de InfectRisk. Assuma que o modelo de 1a. ordem é adquado para modelar essas relações. Obtenha os modelos de regressão ajustados.
  • As funções de regressão estimadas são semelhantes para as quatro regiões? Discuta.
  • Calcule o MSE para cada região. A variabilidade em torno da reta de regressão ajustada é semelhante para as quatro regiões?
  • Construa intervalos de confiança 95% para o coeficiente angular da reta de regressão para cada região. As retas de regressão para diferentes regiões parecem ter mesma inclinação? O que se pode concluir?
  • Construa intervalos de confiança para a resposta esperada correspondendo a InfectRisk = 5, para cada região. O que se pode concluir?
  • Construa intervalos de previsão para um novo hospital em cada região que tenha InfectRisk = 5. O que se pode concluir?

Diagnóstico

  1. Para cada um dos três modelos de regressão ajustados no item (1) da seção anterior, realize o diagnóstico através da análise dos resíduos e apresente um resumo de suas conclusões. O modelo de regressão linear simples clássico é adequado a alguma das situações investigadas?

  2. Ajuste um modelo de regressão linear simples para a variável LengthStay como função de InfectRisk após excluir as observações 47 (X = 6.5 e Y = 19.56) e 112 (X = 5.9 e Y = 17.94). Obtenha intervalos de previsão de 95% de confiança para novas observações que apresentam valores da variável explicativa iguais a 6.5 e 5.9. As observações eliminadas encontram-se nos limites dos intervalos de previsão obtidos? Discuta o significado dos resultados obtidos.

  3. Para os modelos considerando cada região geográfica separadamente no item (2) da seção anterior, realize o diagnóstico através da análise dos resíduos. Todas as regiões aparentam ter mesma variância dos erros? Que conclusões é possível obter a partir da análise?